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Abstract 

Vocal fold (VF) motion is a fundamental process in voice production, and is also a 
challenging problem for direct numerical computation because the VF dynamics depend 
on nonlinear coupling of air flow with the response of elastic channels (VF), which un- 
dergo opening and closing, and induce internal flow separation. A traditional modeling 
approach makes use of steady flow approximation or Bernoulli's law which is known to be 
invalid during VF opening. We present a new hydrodynamic semi-continuum system for 
VF motion. The airflow is modeled by a quasi-one dimensional continuum aerodynamic 
system, and the VF by a classical lumped two mass system. The reduced flow system 
contains the Bernoulli's law as a special case, and is derivable from the two dimensional 
compressible Navier-Stokes equations. Since we do not make steady flow approximation, 
we are able to capture transients and rapid changes of solutions, e.g. the double pressure 
peaks at opening and closing stages of VF motion consistent with experimental data. We 
demonstrate numerically that our system is robust, and models in-vivo VF oscillation more 
physically. It is also much simpler than a full two-dimensional Navier-Stokes system. 



PACS numbers: 43.70Bk, 43.28Ra, 43.28Py, 43.40Ga. 

1 INTRODUCTION 

Vocal folds (VF) are the source of the human voice, and their motion is a fundamental process 

in speech production. In recent years, mathematical modeling of vocal folds has been pursued as 

a viable alternative to direct experimental studies using strobovideolaryngoscopy or electroglot- 

tography techniques. Numerical simulations of VF models then provide us with a valuable tool 
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to understand, monitor and predict various behaviors of normal and disordered voices in vivo. 
Together with models of vocal tract, one can construct a voice simulator which clinicians, speech 
therapists, voice teachers, and otolaryngologists can use to help with their skill improvement, 
diagnosis and patient treatment. 

Since VF motion is mechanical and results from the nonlinear interaction of airflow and 
elastic response of VF, partial differential equations (PDEs) can be written down from classical 
continuum mechanics based on our knowledge of VF structures and air flow characteristics. 
However, the complexity involved is daunting, both in terms of airflow and VF structure, for 
a direct simulation of a complete set of governing equations. Also it is not necessary that one 
needs all the details of such a solution to describe the main VF properties. Modeling effort 
is required to build a smaller set of equations that can capture the essential features of VF 
dynamics. In the past decade, much progress has been made in modeling the elastic aspect 
of VF. There are by now a hierarchy of elastic models for VF, from the two mass model of 
Ishizaka and Flanagan Bogaert 0, to 16 mass as well as the continuum model of Titze and 
coworkers 0, |^, ^ 0. However, the modeling of airflow or the fluid aspect of VF is much less 
explored. 

There are broadly two types of approaches in treating the glottal flow. One approach is to 
combine the Bernoulli's law in the bulk of the flow (steady flow approximation) with empirical 
formulas in boundary layer, flow separation and wake [jl], 0, ^, ^. This approach oversimplifies 
the flow in the sense that PDEs are approximated by algebraic equations. Though 
the approach is a working method for building simple models, it clearly introduces drastic 
approximations. For example, it was realized and concluded that Bernoulli's law is not 



valid during one-fifth of the VF vibration cycle, especially at VF opening and closure. This lack 
of accuracy as a result of deviating significantly from the original PDEs is a major drawback of 
the empirical approach. The other approach is to directly simulate the two dimensional Navier- 
Stokes (NS) system. Two dimensionality is a common assumption for vocal flows p. Alipour et 
al. formulated a steady state simulation with a given glottal geometry |Tl|] . Both ||, |Tl| appeared 
to have been done for fixed channel shape, or in other words, in-vitro VF. To accurately model 
the pulsating nature of the flow during VF vibration however, a time dependent solution is 
more appropriate. Yet in such a case, it is highly difficult to resolve the flows in the presence 
of moving boundaries, closures, and flow separation. Existing works are few in this direction 
although a lot of measurements on the flow characteristics such as intraglottal air pressure and 
flow velocity have been made, by Titze and others [pi llOl O, O, [13, |T5|, [11], ITT . 
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The current status of flow modeling calls for a systematic study of reduced PDE flow models 
and their coupling with existing elastic models. In this paper, we develop an intermediate in- 
vivo PDE model system so that original airflow PDEs are approximated by reduced 
PDEs, not algebraic equations. Consider for now that the two sides of VF are symmetric 
to each other across the centerline, the methodology could be extended to the asymmetric case. 
The air flow is modeled by a quasi-one dimensional (vertical or upward direction) system of 
flow equations. The flow variables (pressure, velocity, density) are averaged quantities over 
the channel cross section of the corresponding ones in two dimensional NS system. Assuming 
that the flow is predominantly in the vertical direction, which is reasonable before flows become 
turbulent in the exit region, we derive the model flow system from the two dimensional isentropic 
compressible Navier-Stokes equations, see the appendix, also [|l^. If the channel is not changing 
in time, the system reduces to the familiar quasi-one dimensional gas dynamic equations in 
studying duct flows in aerodynamics [[191 , POj plj If the channel varies in time, there is 



an additional source term in the flow momentum equation, which turns out to be essential for 
drawing connection with the Titze theory of small VF oscillation pS). One can regard the 



reduced airflow system as a coarse-grained NS system which contains the Bernoulli's law as a 
special case, and inherits time dependent convection mechanisms from the full two dimensional 
NS system. The advantage is that the system is able to handle time dependent effects, such as 
rapid pressure and velocity changes, during VF opening and closing; moreover, it is a lot simpler 
to simulate numerically because all the unknown dependent variables are one dimensional in 
space. Such a system will be coupled to an improved two mass model for the VF cross section 
area to form a complete VF model. The VF cross section area appear as variable coefficients 
in the quasi-one dimensional air flow system. The VF motion is described by how VF cross 
section area varies in time. 

The rest of the paper is organized as follows. In section 2, we introduce the equations of 
the model, and address related modeling issues. In section 3, we discuss numerical method and 
numerical results of model simulation. We show numerically that our model is able to generate 
VF motion in-vivo, and recover several known VF characteristics supported by experimental 
measurements, for example, unequal double pressure peaks at VF opening and closure. We 
also show the robustness of our model by varying subglottal pressure and plotting how air 
volume velocity (air flux) changes as a function of time. Our results reach complete qualitative 
agreement with existing VF flow data. The conclusion is in section 4, and acknowledgement in 
section 5. Section 6 is the appendix on the derivation of our reduced flow model from the two 
dimensional compressible Navier-Stokes system. Table 1, figure captions and figures follow the 
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references. 



2 THE SEMI-CONTINUUM MODEL 

Suppose the larynx is a two dimensional channel with a finite mass elastic wall of cross section 
width A{x,t), and length 2L. The vocal fold is lumped into a sum of two masses connected 
by a spring, and each mass is connected to solid wall by a spring and a damper, the classical 
scenario in the two mass model 0. The air fiows from x = —L to x = L, and is modeled by 



the quasi-one dimensional system derived in ||T8 
• conservation of mass: 



{Ap)t + {puA), = 0, {2.1] 



p air density, u air velocity; 
• reduced momentum equation: 



ipuA)t + (pu'^A), = + A,p + puAt, (2.2) 

p air pressure. 

Assuming that the temperature is maintained as constant, so the airflow is isothermal, then 
the equation of state is: 

p = a'^p, (2.3) 

where a is the speed of sound. The cross section width A is a piecewise linear function in x 
determined by the locations of the two masses (t/i, 1/2), in the classical two-mass model system 
(Bogaert 0, Ishizaka and Flanagan |I|): 

miy'l + riy[ + ki{yi-yo^i) + ki2{yi-y2 + yo,i2) = Fi, (2.4) 
m2?/2 + r22/2 + ^2(2/2 - 2/0,2) + fci2 (2/2 - 2/1 - 2/0,12) = 0, (2-5) 

where Fi = Lg J^^ pdx, Lg the transverse (to the flow) dimension of vocal fold; 2/j's are VF 
openings at locations Xj's, —L < xi < L = X2; Xs = X2 if there is no flow separation, and Xs = 
the location of flow separation if it occurs. The rrii, r^, ki, i = 1,2, are mass density, damping 
and elastic spring constants. Mass one (lower mass) is situated near the VF entrance, and mass 
two (upper mass) is located towards the exit of the glottal region. Following Bogaert f^, Xg will 
be estimated by an empirical formula on the degree of divergence of the VF. Our complete VF 
model is the coupled system (|2.1|) - (^75|) . 

To make the paper self-contained, the derivation of ( ^?Tl) -( p^ is included in the appendix. 
The flow variables in the quasi-one dimensional system are averages over the channel cross 
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section of the corresponding ones in two-dimensional flows. The viscous effect in the flow 
produces the term puAt from the no-slip boundary condition of the two-dimensional flows. 
Without this term, the above system is a familiar one used in gas dynamics (see [^,19, 0, 22 



and references therein) for modeling flows through ducts with variable cross section. It is shown 



T8[| that this extra term is critical in transferring energy from airfiow into the VF, as the 
Titze theory predicted. We have ignored the viscous terms in the momentum equation for 
simplicity, they appear to be higher order. 

Two mass model ( p^.4D -( P3| ) is a recent improvement of the original IF72 []T| in that flow 
separation point is not always at the VF exit, instead it depends on the glottal geometry. Flow 
separation basically refers to a change of flow behavior from being attached to the VF cover via 
a viscous boundary layer to a developed free jet with vortical structures and turbulent wake. 
Because of the vortical buildup, pressure near the wall is typically low, and can be approximated 
by setting it to zero (or ambient pressure) as done on mass two in (|2.5|) when there is no vocal 
tract. In converging glottis, there is no flow separation, however in diverging glottis, it occurs 
if the diverging angle is large enough. It is as yet a challenging problem (no simple theoretical 
prediction) to decide for a flow when and where separation occurs. It is expedient for modeling 
purpose to adopt a working hypothesis supported by experiments 0, |[: 

< 1-1 ^ Xs = X2, (2.6) 

y2/yi>l.l^ x, = xi + ^ — ^^i^, = l.lyi. (2.7) 

10(t/2 - yi) 

Notice that the flow separation location is a variable depending on the diverging angle. It is 
worth pointing out that the assumptions we made for deriving the reduced flow model are all 
valid prior to the separation point. We expect to see a deviation after the flow separation point 
between the reduced flow model and the fully two dimensional NS solutions; however, flow 
pressure post separation is not used in ( |2.5| ). Thus our reduced flow model matches perfectly 
with the improved two-mass model [0. 

We also adopt the elastic collision (stopping) criterion in when the two sides of VF 
approach each other and close. When ?/j's are smaller than a critical level yc, then VF is 
considered closed. Following |I|, ^, {mi,ri,ki) {i = 1,2) are adjusted to closure values. In this 
case, the flow equations are solved only over x G [—L,xi], and in (p.4| )- (|2.5| ) the pressure force 
is adjusted to Fi = Lg f^\pdx. Due to constant input pressure po? pressure at Xi builds up. 
The two mass ODE's (ordinary differential equations) are still running even during VF closure, 
and in due time the increased pressure reopens VF. 

The VF model system is posed as an initial boundary value problem on x G [— L, L], with 
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inlet boundary condition (p, u){—L, t) = (po, Uo), and a zero Neumann type boundary condition 
at exit {px,Ux){L,t) = 0. The advantage of such Neumann type boundary conditions is that it 
helps the flow to go out of the computational domain, which is needed for a stable numerical 
method free of numerical boundary artifacts. Our numerical experiments suggest that the above 
treatment works fine. 

The major difference between our model and that of Bogaert is that we do not make 
quasi-steady approximation on the flow variables, instead we integrate time dependent system 
( [2.1D -( PT^ . This turns out to be particularly important for capturing transients near closure 
and reopening stages of VF motion. 

It is helpful to put the system ( ^31) - 
sound. Then: 

-{Ap)t + {pvA)^ = 0, 
a 

-{pvA)t + [pv^A)^ = -{pA)^ + A^p + pvAt/a, (2.8) 
a 

where typically v = u/a ~ 0.1, the Mach number. If we use the convenient cm-g-ms unit, 
a = Sbcm/ms, 1/a is a small parameter. If we ignore the terms with a, we have exactly 
Bernoulli's law for steady flows. However, these seemingly small terms are essential especially 
during opening stage of VF, and should be kept for an accurate time-dependent solution. 



2.2) into a rescaled form. Let v = u/a, a the speed of 



3 NUMERICAL METHOD AND SIMULATION RE- 
SULTS 

For given VF shape, A{x,t), the flow system (|2.1|) - (|2.2| ) is of the form: 

Ut + {F{U)), = G{U), (3.1) 



so called conservation law (see pj] and references therein) with lower order source term G. The 
function F is the flux function. We implemented a first order finite difference method, where 
time marching is split into two steps. In the first step (t = nk — > (ra + |)fc), we solve the 
conservation law Ut + {F{U))x = with explicit Lax-Friedrichs method [p^ : 

= \{u-_, + f/jVi) - ^ (F(t/;vi) - mU)) ' (3-2) 

where k and h are time step and spatial grid size. Here k must be small enough to ensure 
stability of the difference scheme and to keep the computed flow velocity positive (no back flow 
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is allowed). In step two (t+ | t + fc), we update the solution from U^^2 to f/"+i by implicitly 
integrating ODEs: Ut = G{U) in flow equations, and the two mass equations (2.4)-( p?5D ; where 



we apply central differencing in space and backward differencing in time. In the first step, U 
is updated using VF shape A at time t = nk; in the second step, the ODEs from two mass 
system and source terms are solved to update solutions to {n + l)k. We point out that when 
VF approach closure, the ODE's in step two become rather stiff, and this is the main reason 
to use implicit backward differencing in time [^]. The time step is a variable. It is smaller 
when VFs are closing and the equations are stiff, and is larger when VFs are opening up. 

The parameters used in our calculation are: space grid size h = 0.01125, variable time 
step k G (10^^,10^'^). The time unit is ms, length unit cm, 2L = 0.225 cm, speed of sound 
a = 35cm/ms, uq = Acm/ms, po = 7840 dynes cm^'^. Other two mass model parameters are 
in Table 1. In runs not shown, we have reduced h to half or even smaller sizes, and observed 
similar findings as reported below. 

Now we describe our numerical results, and compare with figures in the literature or those 
from experimental measurements. In Figure |I|, we show a comparison of a cycle of VF vibration. 



The left column is the figure on page 113 of Sataloff's Scientific American article p6[, the right 
column is a plot of our numerically simulated VF vibration cycle. The resemblance is clear. A 
web animation is also available 



In Figure 0, we show our computed air volume velocity (air flux) at the exit of VF, which 
compares well with Fig 6, Fig 7, Fig 8 of Story and Titze ||^. Notice that the airfiux goes 
down to zero gradually at VF closing then drops down abruptly to zero. The drop is due to the 
cutoff introduced in the two mass model for closure, the yc- Below yc, the cross section area A 
becomes very small, and the flow calculation becomes rather stiff. In other runs (not shown), 
we have observed that increasing yc will shorten the curved portion and straighten up the plot 
near closures {t ~ 6, 21, 36 ms). The air volume velocity shows asymmetry, steeper on the right 



side than on the left side, consistent with Fig. 3b of Titze p3| 



Figure ^ is the experimentally measured intraglottal pressure on an excised canine larynx 
from |jl2|] (Fig. 8 there) and [|1^ , which showed the double peak (intraglottal) pressure structure 



respectively at VF opening and closing. Figure |^ is our computed pressure at the grid point 
before lower mass. The double peaks are present and resemble well those in Figure ^ only that 
ours are steeper and higher. Several factors contribute to the difference: (1) we used inviscid 
flow model while there was physical air viscosity in experiments that smooth the solutions, (2) 
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the closure treatment of two mass model differs from the actual VF closure, (3) Figure ^ plots 
the pointwise pressure, not an average pressure over glottis. At the qualitative level however, 
our model solutions are in full agreement with the experimental finding. Notice that the second 
peak is lower than the first. 

To the best of our knowledge, the experimentally observed unequal double pressure peaks 
|T2|, [T^, have not been computed previously in a VF model without coupling to vocal tract. The 
experiment had no vocal tract load. In Story and Titze a computed two peak intraglottal 
pressure plot was given (see Fig 11 [^) using their three mass VF model; however, there is 
additional coupling to vocal tract or an additional subglottal system. Also their computed 
second peak appeared higher than the first peak. The fact that our model can recover the 
experimental double pressure peaks renders strong support for its validity and value. 

We also tested our model robustness under pressure variation. In Figure |^, we show a plot 
of air volume velocity vs time at VF exit for three subglottal pressures: 1584 Pa, 1984 Pa, 2384 
Pa with other parameters same as in Table 1. We see that as subglottal pressures increase with 
other parameters fixed, air volume velocity curves get higher (at peaks) and steeper (at two 
sides). This agrees very well with Fig. 2.14(a), page 78, of K. Stevens j{F^, and is another 
strong support for our model. 

In Figure ^, we show an air particle velocity at VF exit (after upper mass) over three 
vibration cycles for subglottal pressure 2384 Pa, which agrees qualitatively with Fig 3c, page 
1538, of Titze ||2^. For Figure ^ a small amount of additional numerical diffusion is added in 

(O-dD- 



4 CONCLUDING REMARKS 

In this paper, we introduced a new semi-continuum VF model consisting of a reduced PDE flow 
system |jl8| and a recent two mass elastic system 0. The flow part of the model is more physical 



than a traditional treatment with Bernoulli's law yet much simpler than a full two dimensional 
Navier-Stokes system. The reduced PDEs are derivable from the two dimensional compressible 
Navier-Stokes system, and are much more economical for computation. We demonstrated nu- 
merically that the model solutions are in qualitative agreement with known VF experimental 
measurements. In future work, we plan to couple the flow model with more physiological elastic 
VF models, such as H]; compute with higher order flnite difference methods for attaining 
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more accuracy; analyze qualitative properties (phonation thresholds) of model solutions using 
bifurcation methods; and incorporate additional viscous effects in the flow. 
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6 APPENDIX: DERIVATION OF REDUCED FLOW 
SYSTEM 

We derive the fluid part of the model system assuming that the fold varies in space and time 
as A = A{x,t). Consider a two dimensional slightly viscous subsonic air flow in a channel with 
spatially temporally varying cross section in two space dimensions, = ^o(t) = {{x,y) : x G 
[—L, L],y G [— A(a;, t)/2, A{x, t)/2]}, where A{x, t) denotes the channel width with a slight abuse 
of notation, or cross sectional area since the third dimension is uniform. The two dimensional 
Navier-Stokes equations in differential form are (Batchelor page 147): 

• conservation of mass: 

Pt + V ■{pu) = 0; (6.1) 

• conservation of momentum: 

{pu)t = -V- (p(n(g)n)) +div{a-n); (6.2) 
where a is the stress tensor, a = (aij) = —p6ij + dij, and: 

dij = 2p {dj - ^^^y^5jj), Cij = ^{ui^xj + Uj^x,), ixi,X2) = ix,y); 

p is the fluid viscosity; Q{t) is any volume element of the form: {u = (^1,^2)) 

n{t) = {{x,y):xe [a, b] C L],y e [-A{x, t)/2, A{x, t)/2].}. (6.3) 

The equation of state is either polytropic or isothermal. 

The boundary conditions on (p, u) are: 
(1) on the upper and lower boundaries y = ±A{x, t)/2, py = 0, and u = (0, ±At/2), the velocity 
no slip boundary condition; 
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(2) at the inlet, x = —L, p = pq, given subglottal pressure, {ui,U2) = (^^i.o? ""2,0)) given input 
flow velocity. At the exit. {p,Ui,U2)x = 0, to help the waves go out of the domain freely. 

We are only concerned with flows that are symmetric in the vertical. For positive but small 
viscosity, the flows are laminar in the interior of Qq and form viscous boundary layers near the 
upper and lower edges. The vertically averaged flow quantities are expected to be much less 
influenced by the boundary layer behavior as long as A{x,t) is much larger than 0(/i^/^). We 
also ignore effects of possible flow seperation inside Qq when it becomes divergent with large 
enough opening. 

Let us assume that the flow variables obey: 

\ui,y\ ^ \ui,x\, \u2,y\ ^ \ui,x\, away from boundaries of Qq, 
\uy\ ^ \ux\, near the boundaries of Qq, 

\py\ < \px\, throughout Qq. (6.4) 



These are consistent with physical observations in the viscous boundary layers (Batchelor 
page 302), namely, there are large vertical velocity gradients, yet small pressure or density 
gradients in the boundary layers. The boundary layers are of width 0{p^^^). Denote by p, Ui, 
the vertical averages of p and Ui. Note that the exterior normal n = {—Ax/2, 1)/(1 + A^/4)^/^ 
if y = A/2, n = {-Ax/2, -1)/(1 + Al/^fl^ if y = -A/2. 
Let a = X, b = X + 6x, 6x 1, t shghtly larger than to. We have: 

^/ pdV = ^f pJ{t)dV=( ptJ{t)dV+f pJtdV, (6.5) 
dt Jn{t) dt Jn{to) Jn{to) Jn{to) 

where J{t) is the Jacobian of volume change from a reference time to to t. Since Q{t) is now a 
thin slice, J{t) = for small 5x, and Jt = At{t)/A{tQ). The second integral in ( |6.5| ) is: 

/ pJtdV = p^A{to)6x = pAt{t)6x. (6.6) 
Jn{to) A[to) 



The first integral is simplified using (|0 



as: 



/ ptJ{t)dV=[ ptdV = -[ pu-ndS. (6.7) 

jQ.(t„) Jn(t) Jdn(t) 



'n{to) Jn{t) Jdn{t) 

We calculate the last integral of ( |6.7D further as follows 

nA/2 nA/2 



/ pu-nds = / {—pui){x,y,t)dy+ / {pui){x + 6x,y,t) dy 
Jan J-A/2 J-A/2 

+ r^^\-{{), At/2)- {-Ax/2,1) dx 

Jx 
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+ r^^%-{0,~At/2)-{-Aj2,~l)dx 

J X 

= imAt^^^ + y (P A)|,=A/2 + y (P A)|,=-A/2 + 0{{6xf) 

^ {-p-uiA)\l:^^^ + ])At5x + 0{{5xf), (6.8) 

where we have used the smallness of Py to approximate p\y=±A/2 by p and /MT by p ■ liT. 
Combining (H) with: 

^ ( pdV= {pA5x)t + 0((5x)2), (6.9) 
dividing by 5x and sending it to zero, we have: 

{j>A)t + (p ■ UiA)^ = 0, 



which is ( p. 11) . 

Next consider i = 1 in the momentum equation, a = x, 6 = x + 5a;. We have similarly with 



— [ puidV= [ {pui)tdV + [ puiJtdV 
dt Jn(t) Jn(t) Jn(to) 



— puiu-ndS+ / ■ fij dS + pui At6x + 0(6x^). (6.10) 

Jdn(t) Jdn(t) 



idn{t) Jdn{t) 
We calculate the integrals of ( |6.10| ) below. 



4- I pui dV = i-puiA)t5x + Oi{5xf) ^ (p ■ uA)t ■ 5x + 0((5x)2). (6.11) 
dt Jn 

Using ui = on the upper and lower boundaries, a similar calculation as ( |6.8|) gives: 

/ puiu-ndS = {p-ulA)\l+^'' + 0{5xp^/^), (6.12) 

where the smallness of Ui^y in the interior and small width of boundary layer 0(p^/^) gives the 
0(/i^/^) for approximating ulhy ui ■ ui- 

Remark 6.1 Notice that for inviscid flows, we would have an additional term pui At 6x, which 
would cancel the third term on the right hand side of ^6.1(\ ). As a result, the At u/A term would 
he absent from the momentum equation (\2.3i ). 
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Let us continue to calculate: 



+5 

Jdn Jx 

= -pA\l+^^ +pA^6x + 0{{6x] 



2^ 



Noticing that: 



It follows that 



dn = 2/i(ui,^ - + U2,y)/3), d^ = 2fi{ui^y + U2,x)- 

4 _ 2fiAt 
~ 3^'"^'"' 3^4"' 

Thus the contribution from the left and right boundaries located at x and x + Sx is: 

E / , ^11^1 = ^^ir^^ = - ^ir'^ (6.13) 

l,r ''^ 

The contribution from the upper and lower boundaries is: 

E/ duUidS = -duA^6x/2\y=A/2- duAJx/2\y=_A/2 

± Jy=±A/2 

= iJ,Sx^0{dyu)\y=±A/2- (6.14) 
± 

Similarly, 

J2 di2n2dS = iJ,6xJ20{dyu)\y=±A/2- (6.15) 

_t Jy=±A/2 _^ 

Since dyu\y=±A/2 = C'(/i^^/^), the viscous flux from the boundary layers are 0(5x/x^/^), much 
larger than the averaged viscous term S x ^{Aii^^)^ = 0{5xjj). We notice that the vertically 
averaged quantities have little dependence on the viscous boundary layers unless A is on the 
order 0{jj^^'^). Hence the quantities from upper and lower edges in (|6.14| ) and (|6.15|) , and that in 
( |6.12| ), should balance themselves. Omitting them altogether, and combining remaining terms 
that involve only u^, p in the bulk, we end up with (after dividing by 6x and sending it to zero): 



4/i 

(p ■ wA)t + (J> ■ ufA)^ = -{pA)^ + A^p + pulAt + -^{AWx)x - 2pAtj3, (6.16) 



which gives ( p.2| ) in the inviscid limit /i ^ 0. 
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Two Mass Model Parameters in cgs Unit. 

mi = 0.17 g 

777.2 = 0.03 g 

X2 — Xi = 0.2 cm 

xi + L = 0.025 cm 

ki,open = 45 kdynes 

ki,ciosed = 180 kdynes 

yo,i ^ cm 

k2,open 8 kdynes 

k2,dosed = 32 kdynes 

yo,2 = 0.0 cm 

ki2 — 25 kdynes 

yo,12 = C777 

yc = 0.001cm 

A{-L,t) = 2 cm 

ri,open = 17.5 dynes/ {cm s) 

ri,ciosed = 192.4 dynes/ {cm s) 

r2,open = 18.6 dynes/ {cm s) 

r2,ciosed = dynes /{cm s). 
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Figure 2: Simulated VF volume velocity (air flux, cm^ /ms ) vs time at exit of VF from model 
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Figure 3: Experimentally measured intraglottal pressure on excised canine larynx, see Fi^ 
on page 426 of Titze ||l2l . 
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Figure 4: Computed pressure (10^ Pa) change in time at the point right before lower mass. 
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Figure 5: The computed air volume velocity for three values of subglottal pressures. 
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Figure 6: The computed air particle velocity vs. time. 
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